System Dynamics Modelling in R

Agent based modelling

As opposed to system dynamics ABM is a bottom-up approach where individuals (or agents) are allowed to interact with each other autonomously, governed by rules and motivations, with the system emerging from the results of these interactions rather than being formalised in the model itself. This makes ABM very suitable to problems where understanding how individuals behaviours affect systems. ABM is a newer approach

Discrete event simulation

Discrete event simulation is a method which focuses on the concept of entities queuing for a series of ordered events, with each event changing the system in some way (e.g., a conveyor belt moving entities through a well-defined system). This technique is ideal for modelling smaller systems where more detail is required (e.g., patients moving through a hospital). In the below example Huang et al. (2019) explore how to optimise resources to undertake the time critical removal of clots from stroke patients. We can see that patients arrive, forming a queue for resources before being passed into another queue for their next resource. Repeated runs of the model produce slightly different results which can then be analysed to obtain the characteristics of the system.

Examples of DES packages include simmer and PathSimR.

System dynamics

System dynamics by comparison is much more abstract, focusing on continuous flows around a, often large, system rather than the details of individuals or events within the system. Unlike DES, system dynamics treats entities as a fluid controlled by continuous in- and outflows. The top-down methodology of SD is ideal for strategic modelling of the kind of large, interrelated systems common in health care and public health where interventions can have far reaching and unpredictable effects and can take long periods of time to embed.

Summary

In general, if modelling a small, relatively conceptually simple problem over a short period of time with an aim to improve processes DES may be preferred. If modelling large, complex systems, over a long time period where relationships between components of the system and policy changes are of primary interest system dynamics may be preferred. If the aim is to understand how individual behaviour affects systems, especially where behaviour differs across individuals, agent based modelling may be preferred.

Stocks and flows

These are the basic building blocks of a system dynamics models. A stock represents the quantity of a given item at a value in time, for instance the number of patients in hospital on a given day, while a flow is a measure of input or output between stocks, for example the number of admissions or discharges to hospital in a given period. Flows can be thought of as a rate such (e.g., the number of infections occurring per year) while stocks are units (e.g., people, money, hospital beds). Stocks can only change through the flows associated with them.

Auxiliary variables

It is likely that additional variables will be required that may not be stocks or flows, often as a way to limit model complexity. For instance, we may want to define a target or threshold which a stock or flow cannot exceed either as a hard cap or as a fraction of another variable (e.g., to allow the inflow of stock to be proportional to the value of the stock). These are referred to as auxiliary variables.

Feedback loops

Feedback loops are a critical element of system dynamics modelling and refer to circular causal links where the value of a stock influences a flow which in turns changes the stock. Feedback loops are common in systems, for instance if a hospital is full (no spare stock of beds), the inflow of patients admitted will decrease. Because the outflow of patients discharged continues the availability of stock will begin to increase.

A stock can influence a flow directly as in the hospital example, but it can also be modified though intermediate auxiliary variables, for instance we may have a target that the hospital should keep at 85% occupancy and if this is breached the outflow of patients through discharge should be increased, allowing the inflow to continue unabated.

deSolve

As discussed above integration is one of the core components of system dynamics modelling. The deSolve package provides a way to calculate ordinary differential equations (ODE) in R. To use the ode function in deSolve we pass variables containing the stocks, times to simulate over, auxiliary variables, and the model equations, as well as an integration method to use. Running a model would therefore look something like:

data.frame(ode(y = stocks,
               times = simulationtime,
               func = model,
               parms = auxiliaries,
               method = "euler"))

The model should be defined in a function which takes three parameters: the current simulation time, a vector of the current stocks, and a vector of the model parameters. These are wrapped in with and as.list functions to making accessing the parameters easier. The model function returns a list where the first elements is a vector containing the net flows in and out of stocks (referred to as derivatives, see the Integration section). The order these are returned in must match the order they are defined in the stock. In this case we also return a number of additional variables to help understand the model and we rename them for ease of use; these are completely optional.

model <- function(time, stocks, auxiliaries){
  
  with(as.list(c(stocks, auxiliaries)),{
    
    # Auxiliaries
    aEffectiveContactRate <- aContactRate * aInfectivity
    # number of infections created by each infectious person
    
    aBeta <- aEffectiveContactRate / aTotalPopulation
    # per capita transmission rate based on effective contacts per person/step 
    aLambda <- aBeta * sInfected
    # infection rate based on beta and number currently infected
    
    # Flows
    fIR <- sSusceptible * aLambda
    # incidence rate based on the number susceptible and lambda
    
    fRR <- sInfected / aDelay
    # recovery rate based on number infected and recovery delay
    
    # Integral of Stocks
    d_sSusceptible_dt <- -fIR # remove the newly infected from stock
    d_sInfected_dt <- fIR - fRR # get the new number currently infected
    d_sRecovered_dt <- fRR # move those no longer infected to recovered
    
    # Return a list of the important model values
    # The order of the integrals of the stock must match the stock dataframe
    return(list(c(d_sSusceptible_dt,
                  d_sInfected_dt,
                  d_sRecovered_dt),
                IR = fIR,
                RR = fRR,
                Lambda = aLambda,
                CE = aEffectiveContactRate,
                RecDelay = aDelay))
    
  }

DiagrammeR

There are a number of diagram creation packages in R including diagrammeR, ggdag, diagram, and nomnoml and tools such as Visio can be used outside of R and of course technically a model can be built without a diagram, it just becomes harder to understand. For this tutorial we’ll use the diagrammeR package as it is flexible and relatively easy to use. Guidance on using diagrammeR can be found here. Several methods are provided to construct diagrams including piping, dataframes, and create a DOT object (a type of code which creates diagrams). For the purposes of this we will use the DOT language because it allows an (arguably) more logical layout and can allow better integration with the model itself. It can take some playing around to make a nice-looking output.

Other packages

A number of other R packages can optionally be used to make model building easier as well as to enhance the use of the model. Packages within the tidyverse (e.g., dplyr, tidyr) can be used to manipulate data, ggplot, plotly, and echarts4r can be used visualise results with the latter two supporting interactive visuals. Finally Shiny can be used to create interactive dashboards allowing users vary model parameters without writing code. Used together this can be immensely powerful in helping practically develop and use models.

There are a number of packages which can be used to run models developed in Vensim, Stella and other software within R or Python. Examples of these are including a separate script. These have the benefit of being able to develop the model without having to writing the full model as code (which can be tricky) whilst still allowing many of the benefits of using a language like R. If the model was developed in Stella the readsdr package can be used (this can even generate diagrams in R with little code) and if using Vensim pysd2r can be used (note this package also requires Python be installed, more features are available in the base Python package pysd though some advanced features of Vensim such as subscripts may not be supported). More details are given in Appendix: Using other SD software with R.

Stocks

The starting value of stocks is defined at the start of the model and then varies based on the effects of the flows. In this case we need:

• The number of people in the population who are susceptible to the disease at the start of the simulation (Susceptible).
• The number of people in the population with the disease at the start of the simulation (Infected)
• The number of people in the population who have already recovered from the disease at the start of the simulation (Removed)

Auxiliaries

The value of auxiliaries can be set without dependency on the value of other variables (in which case it is termed exogenous because it sits outside the model) or can depend on other variables within the model (and is therefore endogenous). In this case we need:

• The number of contacts between people that lead to passing the disease to another person (Effective Contact Rate). In this case to reduce complexity this will be set outside the model.
• The Total Population is a summation of the Susceptible, Infected, and Removed population. In this case to reduce complexity this will be set outside the model.
• The per capita transmission rate (Beta) of disease from infected people to susceptible people is dependent on the number of effective contacts and the total population.
• The infection rate (Lambda) depends on Beta and the number of people currently infected.
• The Delay between people being infected and being removed from the infected stock from the disease. In our simplified model we will assume that someone is infectious the entire time they are infected. In this case to reduce complexity Delay will be set outside the model.

Flows

While stocks and auxiliary variables are defined as numbers, flows are defined as equations. In this case we need:

• The number of people newly infected is the Incidence Rate. This is an outflow from the Susceptible stock and an inflow to the Infected stock.
• The number of people removed is the Recovery Rate. This is an outflow from the Infected stock and an inflow to the Removed stock.

This is clearly a highly simplified model of disease transmission, not accounting for factors such as difference in sub-populations and changes in policy and behaviour over time, but the first step in making a complex model is to make a simple model. We’ll do this later, but it is always worth considering whether it is worth adding complexity to a model. Complexity will make models harder to understand and maintain and may not provide additional insight into the problem at hand and it is often the case that a relatively simple model provides sufficient understanding to make an informed decision.

Model diagram

The next step is to visualisation our model as a diagram. This both helps with communication to others but also allows us to explore the causal loops and make sure our model makes sense. While some specialised software such as Vensim and Stella are able to create both visuals and models at the same time because R is general purpose these steps must be undertaken separately (though we can provide some linkage which we’ll cover later).

A typical system dynamics diagram looks like:

Vensim SIR Diagram. Duggan (2016).

We’ll go over exactly what each bit of the visual shows later but hopefully, with the definitions of the stocks, flows, and auxiliaries above, the diagram should make some sense and at the very least show how variables relate to each other.

We can make something similar in R:

# ᐩ ⁻ ⭯ ⭮ ⧖ can be used to show loop direction if required (most useful on large diagrams) and flows
# Where variables have a space in the name they should be encapsulated in single quotes

SD_flowchart <- ("
digraph feedback {
    graph []
    
    # Stocks
    node [shape = box]
    Susceptible Infected Removed
    
    # Flows
    node [shape = none]
    'Incidence Rate ⧖' 'Recovery Rate ⧖'
    
    # Auxiliaries
    node [shape = none]
    Lambda Beta 'Effective Contact Rate' 'Total Population' Delay
    
    # Stock Arrows
    Susceptible -> 'Incidence Rate ⧖' [arrowhead = 'none', penwidth = 2]
    'Incidence Rate ⧖' -> Infected [arrowhead = 'normal', penwidth = 2 ]
    Infected -> 'Recovery Rate ⧖' [arrowhead = 'none', penwidth = 2]
    'Recovery Rate ⧖' -> Removed [arrowhead = 'normal', penwidth = 2]
    
    # Flow and Auxiliaries Arrows
    Lambda -> 'Incidence Rate ⧖' [label = ' ᐩ R1⭯', color = blue, arrowsize = 0.5]
    Infected -> Lambda [label = ' ᐩ R1⭯', tailport = 'n', headport = 'e', constraint = false, color = blue, arrowsize = 0.5]
    Beta -> Lambda [label = ' ᐩ', color = blue, arrowsize = 0.5]
    'Effective Contact Rate' -> Beta [label = ' ᐩ', color = blue, arrowsize = 0.5]
    'Total Population' -> Beta [label = ' ⁻', color = blue, arrowsize = 0.5]
    Delay -> 'Recovery Rate ⧖' [label = ' ⁻', color = blue, arrowsize = 0.5]
    Susceptible -> 'Incidence Rate ⧖' [label = ' ᐩ B1⭯', tailport = 'se', headport = 'sw', color = blue, arrowsize = 0.5]
    Infected -> 'Recovery Rate ⧖' [label = ' ᐩ B2⭯', tailport = 'se', headport = 'sw', color = blue, arrowsize = 0.5]
    
    { rank = same; Susceptible; 'Incidence Rate ⧖'; Infected; 'Recovery Rate ⧖'; Removed }
}")

# It may be useful to extract the elements in the diagramme so they can be checked against the parameters in the model to help ensure the diagramme and model reflect each other.
# elements <- read.delim(textConnection(SD_flowchart), header=FALSE, sep=":", strip.white=TRUE) %>%
#   filter(grepl("->", V1))
# 
# elements <- unlist(str_split(elements$V1, pattern = "->"))
# elements <- gsub(pattern = "\\[[^][]*]", replacement = "", x = elements)
# elements <- unique(gsub(pattern = "'| |⧖", replacement = "", x = elements))

DiagrammeR::grViz(SD_flowchart, width = "700px", height = "700px")

Model building

Once we have a concept of how the stock and flows of the system work we can start putting together a model. Firstly, we have to define the model parameters and exogenous variables (i.e., variables where the value is set outside the model such as the number initially infected).

# Create time vector to run the model between. For each step the model is updated.
start <- 0
end <- 30
step <- 0.125

simtime <- seq(start, end, by = step)

# Create stocks vector, with initial values. It can be good practice to prefix all stocks with 's'
stocks <- c(sSusceptible = 99999,
            sInfected = 1,
            sRecovered = 0)

# Create auxiliaries vector, with values. It can be good practice to prefix all auxiliaries with 'a'
auxs <- c(aTotalPopulation = sum(stocks),
          aEffectiveContactRate = 6,
          aDelay = 2)

Next we can write the model. This step is the most complex and it can be useful to refer back to the diagram to ensure the model is specified correctly. From our diagram we know that:

• Susceptible is the prior Susceptible minus the current Infection Rate \(S = S - IR\)

• Infected is the prior Infected plus the current Incidence Rate and minus the current Recovery Rate \(I = I + IR - RR\)

• Recovered is the prior Recovered plus the current Recovery Rate \(R = R + RR\)

• Beta is the Effective Contact Rate divided by the Total Population \(Beta = Effective Contact Rate / Total Population\)

• Lambda is Beta multiplied by the number infected \(Lambda = Beta * I\)

• Incidence Rate is the current Susceptible multiplied by Lambda \(IR = S * Lambda\)

• Recovery Rate is the current Infected divided by the Delay \(RR = I / Delay\)

We can write each of these in R, ensuring that we define them in the necessary order. With these equations written we can place them in the with function (this allows an expression, in this case our model, to be applied to the data), defining the variables we want to return and pass the whole thing as a parameter to the ode function from the deSolve package.

# Write callback function (model equations)
model <- function(time, stocks, auxs){
  
  # Pass the existing stocks and auxiliaries
  with(as.list(c(stocks, auxs)), { 

    # Calculate the updated values 
    
    # Auxiliaries ------
    
    ## per capita transmission rate based on effective contacts per person per step 
    aBeta <- aEffectiveContactRate / aTotalPopulation
    
    ## infection rate based on beta and number currently infected
    aLambda <- aBeta * sInfected
    
    # Flows ------
    
    ## incidence rate based on the number susceptible and lambda
    fIR <- sSusceptible * aLambda
    
    ## recovery rate based on number infected and recovery delay
    fRR <- sInfected / aDelay
    
    # Integral of Stocks ------
    
    ## remove the newly infected from stock
    d_sSusceptible_dt <- -fIR
    ## get the new number currently infected
    d_sInfected_dt <- fIR - fRR
    
    ## move those no longer infected to recovered
    d_sRecovered_dt <- fRR
    
    # Return the new results ------
    return(list(c(d_sSusceptible_dt,
                  d_sInfected_dt,
                  d_sRecovered_dt),
                IR = fIR,
                RR = fRR,
                Lambda = aLambda,
                CE = aEffectiveContactRate,
                RecDelay = aDelay))
    
  })
}

# Call Solver, and store results in a data frame
o <- data.frame(ode(y = stocks, times = simtime, func = model, 
                    parms = auxs, method = "euler"))

This results in a dataframe containing the estimates

Visualise the results

Finally, we can plot the results.

plot_stocks <- ggplot(o) +
  geom_line(aes(x = time, y = sSusceptible, color = "Susceptible")) +
  geom_line(aes(x = time, y = sInfected, color = "Infected")) +
  geom_line(aes(x = time, y = sRecovered, color = "Recovered")) +
  xlab("Day") +
  ylab("Stocks") +
  ggtitle("Example 1: Infectious Disease - Plot of Stocks (People) Over Time") +
  theme_light() +
  theme(legend.position="bottom", legend.title = element_blank())

plot_flows <- ggplot(o) +
  geom_line(aes(x = time, y = IR, color = "Infection Rate")) +
  geom_line(aes(x = time, y = RR, color = "Recovery Rate")) +
  xlab("Day") +
  ylab("Flows") +
  ggtitle("Example 1: Infectious Disease - Plot of Flows Over Time") +
  theme_light() +
  theme(legend.position="bottom", legend.title = element_blank())

grid.arrange(plot_stocks, plot_flows,nrow=2, ncol=1)

Integration

Flows can be plotted as a rate of stock movement over time. Integration is the process we use to calculate what the value of this curve is at a given point in time given an equation. This happens for every stock at every time point. Thus, integration is the mathematical heart of system dynamics.

Unfortunately, however, the most exact method for integration (referred to as analytical integration) is often extremely resource intensive to find and therefore we tend to use a slightly less exact but much faster method called numerical integration. This requires a number of parameters to be set, most importantly the step size and the solver method.

model <- function(t, state, parameters){
        with(as.list(c(state, parameters)), {
                dx <- r*x*(1-x/K)
                return(list(dx))
            })
}

parameters <- c(r = 1.5, K = 10)
state <- c(x = 0.1)

time <- seq(from = 0, to = 10, by = 0.5)
data <- ode(y = state, times = time, func = model, parms = parameters, method = "euler") %>% as.data.frame()

plot05 <- ggplot(data, aes(x = time, y = x)) +
  geom_point(stat="identity", color="#1F78B4", fill="#A6CEE3") +
  stat_function(fun = function(t){0.1*10*exp(1.5*t)/(10+0.1*(exp(1.5*t)-1))}, color = "red") +
  theme_minimal() +
  theme(axis.title.y = element_blank()) +
  ggtitle(paste("Analytical vs Numerical Integration with step size of 0.5"))

time <- seq(from = 0, to = 10, by = 0.125)
data <- ode(y = state, times = time, func = model, parms = parameters, method = "euler") %>% as.data.frame()

plot0125 <- ggplot(data, aes(x = time, y = x)) +
  geom_point(stat="identity", color="#1F78B4", fill="#A6CEE3") +
  stat_function(fun = function(t){0.1*10*exp(1.5*t)/(10+0.1*(exp(1.5*t)-1))}, color = "red") +
  theme_minimal() +
  theme(axis.title.y = element_blank()) +
  ggtitle(paste("Analytical vs Numerical Integration with step size of 0.125"))

time <- seq(from = 0, to = 10, by = 0.0625)
data <- ode(y = state, times = time, func = model, parms = parameters, method = "euler") %>% as.data.frame()

plot00625 <- ggplot(data, aes(x = time, y = x)) +
  geom_point(stat="identity", color="#1F78B4", fill="#A6CEE3") +
  stat_function(fun = function(t){0.1*10*exp(1.5*t)/(10+0.1*(exp(1.5*t)-1))}, color = "red") +
  theme_minimal() +
  theme(axis.title.y = element_blank()) +
  ggtitle(paste("Analytical vs Numerical Integration with step size of 0.0625"))

grid.arrange(plot05, plot0125, plot00625)

model <- function(t, state, parameters){
        with(as.list(c(state, parameters)), {
                dx <- r*x*(1-x/K)
                return(list(dx))
            })
}

parameters <- c(r = 1.5, K = 10)
state <- c(x = 0.1)

time <- seq(from = 0, to = 10, by = 0.5)
data <- ode(y = state, times = time, func = model, parms = parameters, method = "euler") %>% as.data.frame()

plot_euler <- ggplot(data, aes(x = time, y = x)) +
  geom_point(stat="identity", color="#1F78B4", fill="#A6CEE3") +
  stat_function(fun = function(t){0.1*10*exp(1.5*t)/(10+0.1*(exp(1.5*t)-1))}, color = "red") +
  theme_minimal() +
  theme(axis.title.y = element_blank()) +
  ggtitle(paste("Analytical vs Numerical Integration with step size of 0.5 using Euler's"))

time <- seq(from = 0, to = 10, by = 0.5)
data <- ode(y = state, times = time, func = model, parms = parameters, method = "lsoda") %>% as.data.frame()

plot_lsoda <- ggplot(data, aes(x = time, y = x)) +
  geom_point(stat="identity", color="#1F78B4", fill="#A6CEE3") +
  stat_function(fun = function(t){0.1*10*exp(1.5*t)/(10+0.1*(exp(1.5*t)-1))}, color = "red") +
  theme_minimal() +
  theme(axis.title.y = element_blank()) +
  ggtitle(paste("Analytical vs Numerical Integration with step size of 0.5 using LSODA"))

grid.arrange(plot_euler, plot_lsoda)

Delays

It is often the case that the process by which something flows between two stocks takes time. For instance, in the example presented in the last document, there was a delay between people becoming infectious and going on to recover; the greater the delay the more time an individual is infectious, the more contacts they have, and the faster the infection spreads. This example demonstrate how important delays are to system dynamics models.

When setting delays, it is very important to factor in the step size. If we are modelling A&E attendances per day and the length of stay in A&E is 2 hours, then our step size must be at most 2/24 (0.083), and ideally be smaller.

Delays are implemented by dividing the flow by the delay (e.g., flow <- infected / delay). Because the model will run for each time step, a delay of 1 equates to not having a delay.

Changes to systems

Left to its own devices a system dynamics model will merrily model away, but in many cases we want to use external data that can’t simply be included as a model itself or want to introduce some change to how stock flows around the system, either at a given point of time or in response to the model meeting certain conditions. There are two methods we can use to change a system:

• Forcing functions which are ideal when we are modelling a continuous external parameters. An example may be that season affects variables such as the level of admissions or the ability of a virus to spread.
• Events which are ideal to add abrupt changes. Events can occur at a specific time or based on when certain conditions within the model. An example may be adding admissions per day from an external data frame.

start = 0
end = 364
step = end / 12

# a dataframe of known estimates (in this case we generate the values using a sine wave, but we could equally use actual measurements)
exogenous_data <- data.frame(time = seq(start, end, step),
                             transmissibility = sin(seq(0, pi, length.out = length(seq(start, end, step)))))

est_transmissibility <- approxfun(exogenous_data)

approxData <- data.frame(time = 0:364) %>%
  mutate(est_transmissibility = est_transmissibility(time))

ggplot() +
  geom_point(data = exogenous_data,
             aes(x = time, y = transmissibility, color="Actual")) +
  geom_line(data = approxData,
            aes(x = time, y = est_transmissibility, color="Estimate")) 

# Create time vector to run the model between. For each step the model is updated.
start <- 0
end <- 364
step <- 0.125

simtime <- seq(start, end, by = step)

# Create stocks vector, with initial values. It can be good practice to prefix all stocks with 's'
stocks <- c(sSusceptible = 99999,
            sInfected = 1,
            sRecovered = 0)

# Create auxiliaries vector, with values. It can be good practice to prefix all auxiliaries with 'a'
auxs <- c(aTotalPopulation = sum(stocks),
          aEffectiveContactRate = 1,
          aDelay = 7)

model <- function(time, stocks, auxs){
  
  # Pass the exisitng stocks and auxiliaries
  with(as.list(c(stocks, auxs)), { 

    # Calculate the updated values 
    
    # Auxiliaries ------
    
    ## seasonal transmissibility modifier
    
    atransmissibility <- est_transmissibility(time)
    
    ## per capita transmission rate based on effective contacts per person per step 
    aBeta <- (aEffectiveContactRate * atransmissibility) / aTotalPopulation
    
    ## infection rate based on beta and number currently infected
    aLambda <- aBeta * sInfected
    
    # Flows ------
    
    ## incidence rate based on the number susceptible and lambda
    fIR <- sSusceptible * aLambda
    
    ## recovery rate based on number infected and recovery delay
    fRR <- sInfected / aDelay
    
    # Integral of Stocks ------
    
    ## remove the newly infected from stock
    d_sSusceptible_dt <- -fIR
    ## get the new number currently infected
    d_sInfected_dt <- fIR - fRR
    
    ## move those no longer infected to recovered
    d_sRecovered_dt <- fRR
    
    # Return the new results ------
    return(list(c(d_sSusceptible_dt,
                  d_sInfected_dt,
                  d_sRecovered_dt),
                transmissibility = atransmissibility,
                ECR = aEffectiveContactRate,
                Beta = aBeta,
                Lambda = aLambda,
                RecDelay = aDelay,
                IR = fIR,
                RR = fRR
                ))
    
  })
}

# Call Solver, and store results in a data frame
o <- data.frame(ode(y = stocks, times = simtime, func = model, 
                    parms = auxs, method = "euler"))

plot_stocks <- ggplot(o) +
  geom_line(aes(x = time, y = sSusceptible, color = "Susceptible")) +
  geom_line(aes(x = time, y = sInfected, color = "Infected")) +
  geom_line(aes(x = time, y = sRecovered, color = "Recovered")) +
  xlab("Day") +
  ylab("Stocks") +
  ggtitle("Example 1: Infectious Disease - Plot of Stocks (People) Over Time") +
  theme_light() +
  theme(legend.position="bottom", legend.title = element_blank())

plot_flows <- ggplot(o) +
  geom_line(aes(x = time, y = IR, color = "Infection Rate")) +
  geom_line(aes(x = time, y = RR, color = "Recovery Rate")) +
  xlab("Day") +
  ylab("Flows") +
  ggtitle("Example 1: Infectious Disease - Plot of Flows Over Time") +
  theme_light() +
  theme(legend.position="bottom", legend.title = element_blank())

plot_transmissibility <- ggplot(o) +
  geom_line(aes(x = time, y = transmissibility, color = "transmissibility")) +
  xlab("Day") +
  ylab("transmissibility Potential %") +
  ggtitle("transmissibility estimate by day") +
  theme_light() +
  theme(legend.position="bottom", legend.title = element_blank())

grid.arrange(plot_stocks, plot_flows, plot_transmissibility, nrow=3, ncol=1)

# Create time vector to run the model between. For each step the model is updated.
start <- 0
end <- 364
step <- 0.125

simtime <- seq(start, end, by = step)

# Create stocks vector, with initial values. It can be good practice to prefix all stocks with 's'
stocks <- c(sSusceptible = 99999,
            sInfected = 1,
            sRecovered = 0)

# Create auxiliaries vector, with values. It can be good practice to prefix all auxiliaries with 'a'
auxs <- c(aTotalPopulation = sum(stocks),
          aEffectiveContactRate = 6,
          aDelay = 2)

model <- function(time, stocks, auxs){
  
  # Pass the exisitng stocks and auxiliaries
  with(as.list(c(stocks, auxs)), { 

    # Calculate the updated values 
    
    # Auxiliaries ------
    
    ## seasonal transmissibility modifier
    
    atransmissibility <- 1 #est_transmissibility(time)
    
    ## per capita transmission rate based on effective contacts per person per step 
    aBeta <- (aEffectiveContactRate * atransmissibility) / aTotalPopulation
    
    ## infection rate based on beta and number currently infected
    aLambda <- aBeta * sInfected
    
    # Flows ------
    
    ## incidence rate based on the number susceptible and lambda
    fIR <- sSusceptible * aLambda
    
    ## recovery rate based on number infected and recovery delay
    fRR <- sInfected / aDelay
    
    # Integral of Stocks ------
    
    ## remove the newly infected from stock
    d_sSusceptible_dt <- -fIR
    ## get the new number currently infected
    d_sInfected_dt <- fIR - fRR
    
    ## move those no longer infected to recovered
    d_sRecovered_dt <- fRR
    
    # Return the new results ------
    return(list(c(d_sSusceptible_dt,
                  d_sInfected_dt,
                  d_sRecovered_dt),
                transmissibility = atransmissibility,
                ECR = aEffectiveContactRate,
                Beta = aBeta,
                Lambda = aLambda,
                RecDelay = aDelay,
                IR = fIR,
                RR = fRR
                ))
    
  })
}

  # Call Solver, and store results in a data frame
o <- data.frame(ode(y = stocks, times = simtime, func = model, 
                    parms = auxs, method = "euler"))

plot_stocks <- ggplot(o) +
  geom_line(aes(x = time, y = sSusceptible, color = "Susceptible")) +
  geom_line(aes(x = time, y = sInfected, color = "Infected")) +
  geom_line(aes(x = time, y = sRecovered, color = "Recovered")) +
  xlab("Day") +
  ylab("Stocks") +
  ggtitle("Example 1: Infectious Disease - Plot of Stocks (People) Over Time") +
  theme_light() +
  theme(legend.position="bottom", legend.title = element_blank())

plot_flows <- ggplot(o) +
  geom_line(aes(x = time, y = IR, color = "Infection Rate")) +
  geom_line(aes(x = time, y = RR, color = "Recovery Rate")) +
  xlab("Day") +
  ylab("Flows") +
  ggtitle("Example 1: Infectious Disease - Plot of Flows Over Time") +
  theme_light() +
  theme(legend.position="bottom", legend.title = element_blank())

plot_transmissibility <- ggplot(o) +
  geom_line(aes(x = time, y = transmissibility, color = "transmissibility")) +
  xlab("Day") +
  ylab("transmissibility Potential %") +
  ggtitle("transmissibility estimate by day") +
  theme_light() +
  theme(legend.position="bottom", legend.title = element_blank())

grid.arrange(plot_stocks, plot_flows, plot_transmissibility, nrow=3, ncol=1)

# Make a dataframe where we add 10-20 attendances per hour
events_attendances <- data.frame(var = "occupancy",
                                time = 1:(24*7),
                                value = sample(10:20, size = 24*7, replace=T),
                                method = "add")

knitr::kable(head(events_attendances, 5), digits = 2)

model <- function(time, stocks, auxs) {
  with(as.list(c(stocks, auxs)), { 
    dOccupancy <- - aDischargeRate * occupancy
    list(dOccupancy)
  })
}

auxs <- c(aDischargeRate = 0.6)
stocks <- c(occupancy = 0)

times <- seq(from = 0, to = 24*7, by = 1/24)
ed_occupancy <- data.frame(ode(func = model, times = times, y = stocks,
                               parms = auxs, method = "impAdams",
                               events = list(data = events_attendances)))

ggplot(ed_occupancy) +
  geom_line(aes(x = time, y = occupancy, color = "occupancy")) +
  xlab("Hour") +
  ylab("Occupancy of ED") +
  ggtitle("Occupancy with patients added by events") +
  theme_light() +
  theme(legend.position = "none")

# Make a function where we add 10-20 attendances per event but, when the occupancy is greater than 15, we halve the new attendances
eventfun <- function(time, stocks, auxs){
  with(as.list(c(stocks, auxs)), { 
    attendances <- sample(10:20, 1)
    if(occupancy > 15) {attendances <- round(attendances / 2, 0)}
    occupancy <- occupancy + attendances
    return(c(occupancy))
  })
}

model <- function(time, stocks, auxs) {
  with(as.list(c(stocks, auxs)), { 
    dOccupancy <- - aDischargeRate * occupancy
    list(dOccupancy)
  })
}

auxs <- c(aDischargeRate = 0.6)
stocks <- c(occupancy = 0)

times <- seq(from = 0, to = 24*7, by = 1/24)
ed_occupancy <- data.frame(ode(func = model, times = times, y = stocks,
                               parms = auxs, method = "impAdams",
                               events = list(func = eventfun,
                                             time = 1:(24*7))))

ggplot(ed_occupancy) +
  geom_line(aes(x = time, y = occupancy, color = "occupancy")) +
  xlab("Hour") +
  ylab("Occupancy of ED") +
  ggtitle("Occupancy with patients added by events") +
  theme_light() +
  theme(legend.position = "none")

## sourced from: https://stat.ethz.ch/pipermail/r-sig-dynamic-models/2013q1/000160.html

## Derivative function
derivs <- function(t, v, parms) list(c(0, -0.5 * v[2]))

## events
event1 <- function(t, y, parms){
  cat("event1 at t=", t, "\n")
  with (as.list(y), {
    v1 <- v1 + 1
    v2 <- 5 * runif(1)
    return(c(v1, v2))
  })
}

event2 <- function(t, y, parms){
  cat("event2 at t=", t, "\n")
  with (as.list(y), {
    v1 <- v1 + 0.5 * v2
    v2 <- 5 * runif(1)
    return(c(v1, v2))
  })
}

etimes1 <- c(1, 4)
etimes2 <- c(1, 3)
allevents <- sort(unique(c(etimes1, etimes2)))

dispatch <- function(t, y, parms) {
  ## important! conflict resolution
  ##   what if 2 events occur at the same time?
  ## here: use 2nd event
  ##   alternatives: add, average, randomize, ...
  ret <- y
  if (t %in% etimes1) ret <- event1(t, y, parms)
  if (t %in% etimes2) ret <- event2(t, y, parms)
  return(ret)
}

out <- ode(func = derivs,
           y = c(v1  = 1, v2 = 2),
           times = seq(0, 10, by=0.1),
           parms = NULL,
           events = list(func = dispatch, time = allevents))

plot(out)
abline(v=etimes1, col="blue", lty="dashed")
abline(v=etimes2, col="red", lty="dotted", lwd=2)

yini <- c(temp = 18, heating_on = 1, eventoccurs = 0)

temp <- function(t,y, parms) {
  dy1 <- ifelse(y[2] == 1, 1.0, -0.5)
  dy2 <- 0
  dy3 <- 0
  list(c(dy1, dy2, dy3))
}


rootfunc <- function(t, y, parms) {
  time_in_event <- 1
  if(t %in% event_times){ time_in_event = 0}
  yroot <- c(y[1] - 18, y[1] - 20, time_in_event)
  return(yroot)
}

eventfunc <- function(t, y, parms) {
  time_in_event <- 1
  if(t %in% event_times){ time_in_event = 0}
  yroot <- c(y[1] - 18, y[1] - 20, time_in_event)
  whichroot <- which(abs(yroot) < 1e-6) # specify tolerance
  if(whichroot == 2) { y[2] <- 0 } else { y[2] <- 1 }
  if(whichroot == 3) { y[3] <- 1 } else { y[3] <- 0 }
  return(y)
}

event_times <- c(1, 5, 20)
times <- seq(from=0, to=20,by=0.1)
out <- lsode(times=times, y=yini, func = temp, parms = NULL, 
             rootfun = rootfunc, events = list(func=eventfunc, root = TRUE))
plot(out, lwd=2)

Vensim

The following code chunk demonstrates how to load and run an existing Vensim model using pysd2r. More details can be found in the package vignette and the Python pysd package documents. Note some features are unsupported when translating to R (as of writing this includes subscripts).

library(pysd2r)
library(ggplot2)
library(tibble)
library(dplyr)
library(reticulate)

# If required install and connect to the pysd package in Python
# Sys.setenv(RETICULATE_PYTHON = "my_env/bin/python")
# py_install("pysd")

# Connect to pysd
py <- pysd_connect()

#path to model file, in this case held within the pysd2r package
target <- system.file("models/vensim", "Population.mdl", package = "pysd2r") 
py <- read_vensim(py, target)
results <- run_model(py)
results_withtime <- results %>% 
  mutate(TIME = `INITIAL TIME` + cumsum(`TIME STEP`))
ggplot(data=results_withtime)+
  geom_point(aes(x = TIME, y = Population),colour="blue")

Stella

The following code chunk demonstrates how to load and run an existing Sella model using readsdr. More details can be found in the package vignette. Note some features are unsupported when translating to R.

library(readsdr)
library(dplyr)
library(tidyr)
library(ggplot2)
library(deSolve)

filepath <- system.file("models/", "SIR.stmx", package = "readsdr")
mdl <- read_xmile(filepath) 

description <- mdl$description

model_summary <- data.frame(n_stocks = length(description$levels),
                            n_variables = length(description$variables),
                            n_consts = length(description$constants))
print(model_summary)

deSolve_components <- mdl$deSolve_components

all.equal(deSolve_components, xmile_to_deSolve(filepath))

simtime <- seq(deSolve_components$sim_params$start,
               deSolve_components$sim_params$stop,
               deSolve_components$sim_params$dt)

output_deSolve <- ode(y = deSolve_components$stocks,
                      times = simtime,
                      func = deSolve_components$func,
                      parms = deSolve_components$consts, 
                      method = "euler")

result_df <- data.frame(output_deSolve)

head(result_df)

tidy_result_df <- result_df %>% 
  select(time, Susceptible, Infected, Recovered) %>%
  pivot_longer(-time, names_to = "Variable") 

ggplot(tidy_result_df, aes(x = time, y = value)) +
  geom_line(aes(group = Variable, colour = Variable)) +
  theme_classic() +
  theme(legend.position = "bottom")